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LONG-TERM  GOALS 

Our  long-range  goal  is  to  develop  optimization  methods:  1)  to  estimate  the  physical  state  of  the 
ocean  in  order  to  understand  the  present  and  future  conditions  and  associated  variability/uncertainty, 
and  2)  to  utilize  such  forecast  information  for  control-decisions  such  as  optimal  drifter  deployment 
strategy.  This  is  being  accomplished  through  the  use  of  data  assimilation  methods  for  ocean  cir¬ 
culation  models  and  the  study  of  extending  the  assimilation  formulation  to  an  optimal  control 
problem. 


OBJECTIVES 

In  this  effort,  we  study  application  of  the  Monte  Carlo  numerical  techniques  to  problems  of  ocean 
data  assimilation  and  optimal  drifter  deployment.  This  report  focuses  on  our  efforts  on  application 
of  the  particle  filter  to  the  inverse  Lagrangian  prediction  problem  relevant  to  drifter  deployment. 


APPROACH 

An  inverse  Lagrangian  prediction  (ILP)  problem  addresses  retrospective  estimation  of  drifter  tra¬ 
jectories  through  a  turbulent  flow  given  their  final  positions.  In  a  typical  ILP  scenario,  the  launch 
location  of  a  single  or  a  set  of  drifters  in  the  past  is  sought  given  the  present  location(s)  of  these 
drifter(s).  Due  to  chaotic  nature  of  the  forward  Lagrangian  problem  and  limitations  in  accuracy 
and  resolution  of  current  and  wind  data,  it  is  usually  difficult  to  expect  an  unique  and  determinis¬ 
tic  answer  to  an  ILP  problem.  It  may,  however,  be  possible  to  estimate  the  launch  site  and  time 
statistically  so  that  the  drifters  deployed  in  the  estimated  region  and  time  would  have  the  largest 
probability  of  arriving  at  the  desired  final  location.  For  most  practical  problems  it  is  desirable  to 
minimize  the  optimal  deployment  region  while  maximizing  the  probability  of  successful  delivery. 

One  approach  to  solving  the  ILP  problem  is  to  simulate  an  ensemble  of  Lagrangian  trajectories 
backwards  in  time  using  the  known  final  locations  and  a  stochastic  model  of  the  flow  field.  Due  to 
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the  typically  fast  rate  (exponential  or  geometrical  functions  of  time)  dispersion  of  the  trajectories, 
however,  the  distribution  of  the  drifter  locations  tends  to  be  too  diffuse  to  be  able  to  locate  the 
launch  site. 

We  have  investigated  a  numerical  method  that  employs  the  particle  filter  to  control  the  spread  of 
the  drifter  location  (Chin  and  Mariano,  2008).  The  particle  filter  is  a  general  sequential  estimation 
method,  is  highly  flexible,  and  is  based  on  Monte  Carlo  simulations  of  the  state  trajectory  (Doucet 
et  al,  2001).  Application  of  the  particle  filter  to  general  data  assimilation  problems  has  been  under 
consideration  (van  Leeuwen,  2003).  Only  a  special  case,  known  as  the  ensemble  Kalman  filter 
(Evensen,  1994)  which  is  based  on  formulas  optimal  only  under  the  assumption  that  the  state 
(prognostic)  variables  have  a  jointly  Gaussian  distribution,  has  so  far  been  found  to  be  practical, 
because  the  particle  filter  otherwise  demands  a  prohibitive  number  of  samples  of  ocean  state  to 
be  effective.  However,  since  simulation  of  Lagrangian  trajectories  can  be  repeated  easily  without 
a  prohibitive  demand  on  computational  resources,  an  unrestricted  particle  filter  (that  allows  non- 
Gaussian  distributions)  would  be  applicable  to  ILP  problem. 

For  ILP,  we  assume  knowledge  of  an  empirical  time  characteristics  of  ensemble  spread ,  quantified 
here  by  the  standard  deviation 
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where  r(f)  is  the  unknown  drifter  trajectory,  r n(t)  is  the  nth 

sample  of  such  trajectory  by  sim¬ 
ulation,  and  r (t)  =  J2n=irn(t)/N  is  the  ensemble-mean  of  such  samples.  In  ILP,  the  drifter 
trajectories  simulated  backward  in  time  are  expected  to  converge  towards  each  other  due  to  causal¬ 
ity.  For  example,  in  our  test  cases,  we  have  derived  Dr(t)  oc  tc  for  c  ~  4/3  by  forward  simulations. 
To  formally  express  the  information  with  which  to  constrain  the  backward  trajectory  ensemble,  we 
let  s (t)  be  a  fictitious  “noisy  observation”  of  the  unknown  trajectory  r(t) 


s  (t)  =  r  (t)  +  eft)  (2) 

where  off  )  is  vector  of  random  observation  errors  each  with  a  known  variance  E2.  Assuming  that 
r(f)  and  eft)  are  uncorrelated,  the  variance  of  s (t)  would  become  Dr(t )2  +  E2.  Since  the  mean 
of  s (t),  or  an  observation  of  the  mean  trajectory,  is  not  available,  we  estimate  it  in  a  bootstrapping 
fashion  using  the  ensemble  mean  r(t)  of  the  on-going  simulation.  The  probability  density  function 
(PDF)  ps|r  of  the  observation  s  conditioned  on  the  unknown  r  is  used  by  the  particle  filter  algorithm 
to  constrain  the  state  trajectory;  the  specific  algorithm  is  called  the  resampled  particle  filter  (RPF; 
Chin  et  al  2007).  By  experimentation,  we  have  found  that  an  alternative,  non-Gaussian  PDF  is 
effective  in  imposing  the  constraint  onto  the  ensemble  of  backward  trajectories: 
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where  c  is  a  normalization  constant  and  F  is  a  constant  parameter  to  control  “flatness”  of  the  PDF. 
For  F  —  1,  ps|r  would  become  a  Gaussian  PDF.  We  use  F  —  3  so  that  the  PDF  would  have  a 
relatively  flat  peak  near  its  maximum.  Choosing  the  larger  value  of  F  would  give  more  equal 
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importance  (probability)  to  the  ensemble  members  within  a  certain  distance  from  the  maximum, 
rather  than  favoring  those  in  the  immediate  vicinity  of  the  maximum.  More  complete  details  of  the 
numerical  procedure  can  be  found  in  (Chin  and  Mariano,  2008). 


WORK  COMPLETED 

1)  Performance  of  particle  filter-based  ILP  solution  method  has  been  evaluated  using  test  experi¬ 
ments  of  “source  estimation”  and  “array  deployment”  (see  Result  below). 

2)  Inter-comparison  study  involving  the  EnROIF  assimilation  system,  a  Monte-Carlo  (or  ensemble- 
based)  enhancements  to  the  existing  ROIF  method  (Chin  et  al  2002),  is  complete  for  the  1/12- 
degree  resolution  Gulf  of  Mexico  circulation  model  using  HYCOM,  and  a  report  (Srinivasan  et  al, 
2008)  is  being  finalized. 


RESULTS 

To  examine  performance  of  the  particle  filter-based  ILP  solution  method,  two  types  of  controlled 
experiments  have  been  conducted.  In  source  estimation  experiment,  the  final  locations,  called 
targets  and  denoted  as  Xm,  are  given  by  forward  trajectory  simulations  using  a  single  launch 
site  (for  five  different  cases  shown  in  Fig.  1).  The  goal  is  to  estimate  the  launch- site  probability 
distribution  given  the  target  locations  using  the  filtered  backward  trajectory  ensemble.  In  array 
deployment  experiment,  the  target  locations  are  arbitrarily  chosen  and  there  is  no  guarantee  that  a 
single  launch  site  is  sufficient  to  reliably  deliver  drifters  to  all  targets. 

We  evaluate  benefit  of  the  particle  filter  (or  more  specifically  RPF)  by  comparing  two  ensembles  of 
trajectories.  One  is  an  ensemble  produced  with  the  RPF  procedure  denoted  as  rPPF,  n  —  1, . . . ,  N; 
the  other  is  an  ensemble  of  trajectories  without  any  constraint  and  denoted  as  rFns.  To  compare  the 
two  ensembles,  the  launch  site  distribution  estimated  by  each  ensemble  is  used  to  initialize  some 
test  drifters  for  forward  trajectory  simulations.  The  target  locations  estimated  by  the  test  drifters 
can  then  be  used  to  evaluate  statistical  accuracy  in  reproducing  the  known  target  locations  Xm, 
m  —  1, . . . ,  M. 

Two  skill  scores  are  computed  to  compare  the  ensembles  rPPF  and  rFns.  By  letting  G'(Xm)  be  the 
chance  (in  %)  of  the  mth  target  being  reached  by  any  of  test  drifters,  we  define  the  coverage  score 
to  be  7  =  minm  G(Xm).  We  also  define  the  //  to  be  average  chance  (in  %)  of  a  test  drifter  to  reach 
any  of  the  targets.  A  higher  p  value  indicates  that  a  drifter  from  the  ensemble  is  less  likely  to  miss 
a  target  and  that  the  drifter  destination  is  more  likely  to  be  focused  near  a  target  location.  We  hence 
call  //  the  resolution  score. 

For  the  source  estimation  experiment,  the  area  bounded  by  the  95-percentile  density  contours  (thick 
lines,  Fig.  2)  of  rFns  expands  as  the  simulation  progresses  as  expected,  while  the  analogous  area  for 
rPPF  contracts  in  accordance  with  Dr(t).  By  the  launch  time,  rPPf  has  surrounded  the  launch  site 
(star,  lower-right  panel)  tightly  with  the  10-percentile  contour,  while  the  rFns  could  surround  the 
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launch  site  only  loosely  with  the  50-percentile  contour  (lower-left  panel).  This  observation,  that 
the  launch  site  estimate  is  significantly  tighter  for  the  filtered  ensemble  rPPF  than  unconstrained 
ensemble  rFns,  is  generally  true  for  each  of  five  launches  that  we  tested  (Fig.  3).  The  area  /1PPF 
enclosed  by  the  95-percentile  contours  of  rPPF  was  19  ~  32%  of  the  unconstrained  counterpart 
AF"S  (Table  1).  The  coverage  scores  were  perfect  (7  =  100)  for  both  rFns  and  rRPF  in  all  five  cases, 
implying  that  the  filtered  solutions  are  just  as  effective  in  delivering  drifters  to  all  target  sites.  The 
resolution  scores  are  1.6  to  1.9  times  higher  for  the  filtered  results  than  their  unconstrained  coun¬ 
terparts,  indicating  that  1.6  to  1.9  times  more  drifters  launched  according  to  the  filtered  solution 
would  reach  a  target  than  those  prescribed  by  the  unconstrained  solution. 

In  the  array  deployment  experiment,  the  launch  site  for  each  of  four  arrays  of  5  x  5  targets  (Fig.  4) 
is  estimated.  Again,  the  filtered  ensemble  rPPF  has  yielded  more  compact  estimate  of  the  launch 
site  than  the  unconstrained  ensemble  rFns,  as  indicated  by  the  ratio  /1PF’F//1F5 s  ranging  11  ~  23% 
in  four  cases  (Table  2).  However,  while  the  unconstrained  solution  has  perfect  coverage  scores, 
the  filtered  solution  has  missed  perfect  coverage  in  all  but  one  cases  (7RPF,  Table  2).  In  particular, 
there  were  one  target  site  each  in  the  “NW”  and  “SC”  array  cases  and  two  target  sites  in  the  “SW” 
case  that  had  less  than  perfect  (G(Xm)  <  100)  coverage.  To  remedy  the  coverage  issue,  we  have 
performed  the  supplemented  RPF  (SPF)  procedure,  where  extra  simulations  of  backward  drifters 
are  initialized  only  at  the  target  site(s)  Xm  with  imperfect  coverage  G(Xm)  <  100.  The  resulting, 
enhanced  solution  still  is  nearly  as  compact  (Fig.  5)  as  before  (A|PF/A|51S,  Table  2),  yet  now  with 
perfect  coverage  (7SPF,  Table  2). 


IMPACT/APPLICATIONS 

We  have  explored  a  particle  filter  approach  to  solve  the  inverse  Lagrangian  prediction  problem 
by  an  ensemble  simulation  of  backward  trajectories.  The  numerical  experiments  demonstrate  that 
ensemble  spread  can  be  controlled  using  a  constraint  derived  empirically  and  that  the  constrained 
solution  leads  to  a  spatially  more  compact  estimate  of  the  launch  site.  The  constrained  solution 
is  thus  more  efficient  than  the  unconstrained  counterpart,  while  not  compromising  much  effec¬ 
tiveness  in  delivery  to  the  intended  target  sites.  Due  to  high  demands  for  shipping  resources  in 
drifter  deployment,  adopting  the  technique  to  actual  operations  and  evaluating  its  benefits  would 
be  potential  topics  of  future  investigation. 

The  particle  filter,  more  specifically  the  resampled  particle  filter,  is  well  suited  for  realization  of 
the  constrained  trajectory  simulation  due  to  flexibility  of  the  method.  While  the  particle  filter 
methods  are  known  generally  to  require  a  large  number  of  samples  to  perform  well  (Chin  et  al 
2007),  this  issue  should  not  become  a  limiting  factor  in  Lagrangian  applications  due  to  relatively 
low  computational  cost  of  simulating  a  trajectory  given  a  background  flow  field. 

Systematic  errors  in  the  background  flow  would  pose  the  main  challenge  in  application  of  the 
method  studied.  The  main  issue  is  not  that  these  errors  can  affect  the  mean  trajectory  of  a  drifter 
ensemble  (and  hence  the  estimate  of  the  launch  site)  but  that  there  is  general  lack  of  viable  models, 
statistical  or  otherwise,  of  such  errors.  Use  of  basis  functions  such  as  empirical  orthogonal  func¬ 
tions  may  be  among  the  few  available  techniques  to  characterize  systematic  flow  errors.  Efforts  are 
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under  way  to  enhance  the  presented  techniques  to  incorporate  some  assumed  mathematical  forms 
for  systematic  errors  in  the  background  flow  field. 


RELATED  PROJECTS 

The  data  assimilation  (EnROIF)  component  of  this  project  has  been  associated  with  the  U.S.  GO- 
DAE:  Global  Ocean  Prediction  with  the  Hybrid  Coordinate  Ocean  Model  (HYCOM),  in  collabo¬ 
ration  with  the  HYCOM  Consortium  (http://hycom.rsmas.miami.edu). 

The  ILP  solution  methodology  developed  in  this  project  may  see  applications  in  larval  dispersal 
study  (Cowen  et  al,  2008)  using  HYCOM  simulated  velocity  fields  over  the  Intra- America  Seas 
region. 
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Table  1.  Skill  scores  from  the  source  estimation  experiment. 
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Far  East 

0.316 

100 

100 

49.4 

80.9 

1.64 

Center  Saddle 

0.240 

100 

100 

33.1 

64.2 

1.94 

Table  2.  Skill  scores  from  the  array  deployment  experiment. 
“ SPF ”  denotes  supplemented  RPF  runs. 
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Figure  1.  Left:  The  five  launch  sites  (stars)  and  the  background  flow  at  the  launch  time  (dark 
contours  are  anti-cyclonic;  light  contours  are  cyclonic)  for  the  source  estimation  experiment. 
The  launch  sites  are  named,  counter-clockwise  from  upper-left,  “NW  Saddle ”,  “ Jet ”,  “Cold 
Core ”,  “Far  East”,  and  “Center  Saddle”.  Right:  The  30  target  sites  (dots)  corresponding  to  the 
“Jet”  launch  site  (star). 
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Figure  2.  Progression  of  the  density  of  backward  drifters  initialized  at  each  of  the  30  target  sites 
shown  in  Fig.  1  ( right panal).  The  95 percentile  (thick  lines)  as  well  as  50  and  10 percentile  (thin 
lines)  contours  of  the  drifter  density  are  shown  after  10,  20,  and  30  days  of  simulation.  The  left 
and  right  columns  show  the  simulations  without  and  with  the  particle  filter,  respectively.  The 
source  (star)  is  located  more  accurately  with  the  particle  filter  after  the  30-day  analysis  period. 
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Figure  3.  The  95  percentile  contours  of  the  drifter  density  after  30  days  of  inverse  simulation, 
estimating  the  four  of  five  launch  sites  shown  in  Figure  1  (see  Figure  2  for  the  one  remaining 
launch  site).  Results  obtained  with  ( dark  contour )  and  without  (light  contour)  filter  constraint 
are  shown  in  each  panel.  Again,  simulations  with  the  particle  filter  can  localize  each  launch  site 
more  tightly. 
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Figure  4.  The  four  5x5  arrays  of  target  sites  and  the  background  flow  at  launch  time  (dark 
contours  are  anti-cyclonic;  light  contours  are  cyclonic)  in  the  array  deployment  experiment. 
The  target  arrays  are  named,  clockwise  from  top-left,  “NW”,  “NE”,  “ SC ”,  and  “SW”. 
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Figure  5.  The  95  percentile  contours  of  drifter  density  indicating  the  potential  launch  site  to 
deploy  the  shown  drifter  array.  Each  panel  shows  results  obtained  with  ( dark  contour)  and 
without  (light  contour)  the  particle  filter. 
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